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ABSTRACT 

Particles ejected from the Sun that stream through the surrounding plasma of the solar wind are causing instabilities. These generate 
wavemodes in a certain frequency range especially within shock regions, where particles are accelerated. The aim of this paper 
is to investigate of amplified Alfvenic wavemodes in driven incompressible magnetohydrodynamic turbulence. Results of different 
heliospheric scenarios from isotropic and anisotropic plasmas, as well as turbulence near the critical balance are shown. The energy 
transport of the amplified wavemode is governed by the mechanisms of diffusion, convection and dissipation of energy in wavenumber 
space. The strength of these effects varies with energy and wavenumber of the mode in question. Two-dimensional energy spectra of 
spherical k-space integration that permit detailed insight into the k\\ /^-development are presented. 

The evolution of energy injected through driving shows a strong energy transfer to perpendicular wavemodes. The main process 
at parallel wavemodes is the dissipation of energy in wavenumber space. The generation of higher harmonics along the h\\ axis is 
observed. We find evidence for a critical balance in our simulations. 

Key words. Magnetohydrodynamics (MHD) - Turbulence, Sun: coronal mass ejections (CMEs) - Sun: particle emission 



1. Introduction 

The Sun emits energetic particles during coronal mass ejections 
(CME) and solar flares. The observed energies connected to 
these particles reach GeVs in extreme cases. At first, solar ener- 
getic particle (SEP) events were linked exclusively to solar flares 
(IForbushj 1946), but once CMEs were first observed in the sev- 
enties, it became clear that they have a major role in the gene- 
sis of SEP events ( |Kahler et al.|197 8') The CME-driven shocks 
are now believed to be the major acceleration agent during the 
strongest SEP events (e.g., Reames 1999) and the most plau- 
sible acceleration mechanism in operation is diffusive shock ac- 
celeration (DS A) ( |Bell T978 ) For further reading into the matter 
observations and models of CMEs we recommend Forbe s~et al.l 
|fl2006|ll and |Chen[p)nf 1 

In DSA, particles are accelerated through repeated cross- 
ings of the shock compression front, each crossing giving a 
small boost to the particle energy. The shock crossings are me- 
diated through interactions of particles with background waves. 
Furthermore, the particles amplify wavemodes. This yields a 
coupled system of waves and particles as described, e.g., in 
Vainio & Laitinen (2007) Since DSA is a resonant wave-particle 
process, it is now interesting to see if energy transfer between 
wavemodes will affect different particle energies. A detailed 
treatement of the DSA process itself is given in |Drury| (fT983|l j, 
and Schlickeiser ( 2002) To enable acceleration to the GeV en- 
ergies through DSA, the upstream medium needs to be very 
turbulent. The ambient solar wind turbulence levels are gen- 
erally too low to account for this acceleration mechanism to 
operate beyond the MeV regime (Vainio 2006), but accelera- 
tion to the highest energies can still proceed through the strong 
amplification Alfven waves in the ambient medium through 
streaming instabilities driven by the accelerated particles them- 
selves. Analytical ( |Lee|2 005) and numerical models of the self- 
consistent particle acceleration in coronal plasma have been pre- 
sented ( |Vainio & Laitin enf2007 Ng & Reames|2 008 ) showing 
that the wave generation process is strong enough to account for 



the turbulence responsible for fast scattering of particles from 
one side of the shock to the other. The models also show that 
if particles are accelerated to hundreds of MeVs, the waves will 
grow to nonlinear amplitudes close to the shock. 

It is crucial to understand the processes that govern the tur- 
bulent waves because the mechanism of DSA is strongly con- 
nected to the wave-particle interactions. Nonlinear Alfven waves 
may interact with each other, which may lead to three impor- 
tant effects from the point of view of particle scattering: firstly, 
wave decay through three-wave interactions may limit the wave 
amplitudes in the shock environment; secondly, the wavenum- 
ber spectrum of the Alfven waves may be altered so that the 
waves fall out of resonance with the particles; thirdly, the cross- 
helicity state of the resonant waves may also change, which af- 
fects the scattering center compression ratio of the shock and 



thus the accelerated particle spectrum (Vainio & Schlickeiser 
1998) According to our knowledge, nonlinear wave transport 
models have not yet been applied to the SEP acceleration prob- 
lem. 

In this paper we will take the first steps towards understand- 
ing the nonlinear evolution of waves generated by energetic par- 
ticle beams in the solar corona. We concentrate on a simpli- 
fied model, where the beam-generated wave component is rep- 
resented by a Gaussian peak in parallel wavenumber that fol- 
lows the interaction of this spectral component with a quasi- 
isotropic background turbulence driven randomly in an incom- 
pressible magnetohydrodynamic (MHD) simulation. The tur- 
bulent plasma environment mimics the fast solar wind, where 
Alfven waves are observed to be the dominant species ( |Telloni| 



et al. 2009; Bruno & Carbone 



can be assumed. It is especial 



2005 i , hence incompressibility 



y interesting to see the energy 



transport in parallel and perpendicular direction. Taking into ac- 
count the resonance condition, only energy transport parallel to 
the background magnetic field will alter the transport of parti- 
cles with an energy different from the incident particle's energy. 
On the other hand, most turbulence theories for incompressible 
plasmas predict perpendicular energy transport. 
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2. Numerical model 

Despite MHD-turbulence has been studied for roughly 70 years 
after it was initiated by Hannes Alfven, it is still a controversial 
topic. We focus on incompressible turbulence for which some 
promising progress has been made in recent years. 

One commonly observed property is the characteristic en- 
ergy spectrum E(k) following a power-law with a slope of -5/3, 
which is commonly referred to as the Kolmogorov-type spec- 
trum. Kolmogorov (1941) predicted this power-law for hydro- 
dynamic turbulence by using dimensional analysis and scaling 
behaviour assuming isotropy. The basic system of turbulence 
evolution has also been explained by Kolmogorov. On large 
scales, i.e. small wavenumbers, energy is injected into the turbu- 
lent fluid. This energy decays by generating smaller structures 
up to the smallest scales where dissipation becomes dominant. 
Consequently, dissipation maintains the energy flow from small 
towards large k. This is also the reason for dividing the spec- 
trum into driving-, inertial-, and dissipation range. Although 
Kolmogorov model of turbulence was first discussed in connec- 
tion with neutral fluids, it seems to be valid in the magnetohy- 
drodynamic case as well. 

A different approach to Kolmogorov's theories by Iros hnikov] 



1964) and Kraichnan (1965 ) assuming a local mean magnetic 



field and assuming Alfven wave packets led to an exponent of 
-3/2. The problem of the Iroshnikov-Kraichnan (IK) model is 
the assumption of isotropy because a background magnetic field 
will lead to a preferential direction in space caused by wave in- 
teraction resonances. The IK model implies resonant three wave 
interactions within a weak turbulent regime. In magnetised in- 



compressible plasmas these interaction rates are empty (Sridhar 
|& Goldreich| ( fT994l )). 

Goldreich and Sridhar (GS) describe anisotropic turbulence 



and distinguish between its weak (Sridhar & Goldreich||1994|l| 
and strong state (Goldr eich & Sridhar|1995| Their assumption 
of strict separation between three and four wave interactions 
was controversially discussed and an intermediate state was in- 
troduced (Goldreich & Sridhar 1997) These theories describe 
Alfvenic turbulence evolution towards the perpendicular direc- 
tion. In the weak turbulence regime four wave interactions are 
the underlying process, referring to the GS-framework. Due to 
the resonance condition energy transfer to parallel wavenum- 
bers is not possible. It is still debated wheter four wave in- 
teractions are indeed the basic mechanism in weak turbulence 
is questionable. In recent theories the intermediate turbulence 
(Goldreich & Sridhar (1997)), which is based on three wave 
interactions, replaces the weak four wave interaction model 
( |Lithwick & Goldre ich (2003)). However, strong turbulence is 
dominated by nonresonant three wave interactions, which leads 
to an anisotropic energy-cascade in the perpendicular direction. 
Parallel evolution is not caused by cascading. One of the main 
achievements of the Goldreich and Sridhar theory is that is ex- 
plains the Kolmogorov-type energy spectrum for an anisotropic 
regime. This could explain the observed -5/3 slope in parts of 
the solar wind ( |Bruno & Ca rbone 2005 ) where Kolmogorov- 
theory is not applicable. 

The region of the heliosphere we are interested in is within 
the weak turbulence regime, with magnetic field fluctuations de- 
fined as 



dB = B - Ba 



(1) 



with a mean value of (dB) = 0, which leads to (B) ~ It is 
observed that the solar wind magnetic fluctuations decrease as 
dB 2 oc r~ 3 , while the background field decreases as oc r~ 4 



( jBavassano et al.|1982||Bruno & Carbone|2005| Consequently, 
the dB/Bo ratio within the heliosphere ratio is increasing with 
distance to the Sun (Hollweg et al.|2010| ) 

A remark on the notation, the magnetic background field Bo 
is defined towards the z-direction within our simulations and 
hence is also noted as B^e z . The parallel direction is, therefore, 
defined as the z-directions and the x- and y-direction are the per- 
pendicular directions. For symmetry reasons there will be no fur- 
ther distinction between the two perpendicular directions and all 
plots show values averaged over the azimuthal angle in cylindri- 
cal coordinates of x and y. 

For small perpendicular wave numbers the transport in per- 
pendicular direction is dominating until k ± v ± is of the same or- 
der as the Alfven cascading time k\\v A - This means that in addi- 
tion to the perpendicular cascade, a cascade in the parallel di- 
rection will occur as well and energy will be transferred towards 
smaller parallel spatial scales. Accordingly, the parameter 



k±v±_ 
k\\v A ' 



(2) 



where v A is the Alfven velocity, is of the order of unity. This 
state is called the critical balance and was first introduced by 
Goldreich & Sridhar| ( |1995| ) In this state the linear wave period 



of the Alfven waves are comparable to the intrisically nonlin 
ear timescale. If £ ~ 1, the fluctuations become more correlated 
along the parallel direction, up to l\\ ~ v A /(k ± v ± ) as indicated by 
Eq. [2] Then the turbulence is clearly within the strong regime. 
This means that the fluctuations become comparable to Bo and 
the nonlinear term is not small anymore (Perez & Boldyrev 
[20081 ) 

High Reynolds numbers in combination with massive energy 
injection, as seen in, e.g., the solar wind, are strong indicators of 
a highly turbulent state. In situ measurements of the energy spec- 
trum ( |Tu et al.| 1989 ) agree with this fact. To simulate conditions 
within the turbulent heliospheric plasma, the research group at 
the University of Wiirzburg has developed a simulation code, 
Gismo. 

Gismo is an incompressible pseudospectral MHD-code that 
is fully parallelised and capable of efficiently using massive 
computing clusters. The basis of the simulation software is to 
solve the following set of incompressible MHD-equations: 



du 

lit 
db 



b-Vb-u-Vu-VP + v,,V m k 



= b-Vu-u-Vb + j]V 2h b 



dt 

V-h =0 
V ■ b = 0, 



(3) 



where b = B/ yjA-ng is the normalised magnetic field, u is the 
fluid velocity and g is the constant mass density. The diffusion 
coefficient related to viscous and Ohmic dissipation is denoted 
by v v and rj. A common approach in pseudospectral methods is 
to amplify the diffusion term by a power of h - resulting in hy- 
perdiffusivity. This artificial enhancement of the dissipation is 
necessary to reach a saturated state of turbulence within a rea- 
sonable timescale. It is a methodic problem, since pseudospec- 
tral approaches do not strongly suffer from dissipative numerical 
effects. The only intrisic energy loss of the system is caused by 
antialiasing, which we discuss below. Furthermore, the parame- 
ter v is introduced as a global diffusivity with rj = v v = v. Hence 
magnetic resistivity and viscous damping are not distinguishable 
anymore. This is the case for the magnetic Prandtl number of the 
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order of unity, which is valid within the regime of Alfven wave 
turbulence where an equipartion between magnetic and kinetic 
energy can be assumed (Bershadskii 2002; Bigot et al. 2008) 
The pressure term VP fulfills the closure condition for incom- 
pressibility ( Maron & Goldreich|2001|) | 



V 2 P = Vb : Vb - Vh : Vh. 



(4) 



These equations are solved in Fourier space by using pseu- 
dospectral methods that lead to the componentwise equations 



at 

dbg 
k a b = 0, 



kfykp I / \ , 

= -iky \8 a p - —j- J [upUy - bpb y j - vk 'u a , 



= -ikp (ujjb a - bpu a ^ - vk 2h b a , 



(5) 



where the tilde-notation stands for quantities in Fourier space. 
The components of the wavevector are written as k a . 

In the incompressible regime of a magnetised plasma the 
MHD-turbulence consists of only two types of waves, which 
propagate along the parallel direction - the so-called pseudo- 
and shear Alfven waves. The Former are the incompressible 
limit of slow magnetosonic waves and play a minor role within 
anisotropic turbulence (Maron & Goldreich (2001 )). The pseudo 
Alfven waves polarisation vector is in the plane spanned by the 
wavevector k and Bo- The shear waves are transversal modes 
with a polarisation vector perpendicular to the k - Bo plane. 
They are circularly polarised for parallel propagating waves. 
Both species exhibit the dispersion relation a> 2 — (v A k\\) 2 . Note 
that the shear mode seems to be dominant because pseudo waves 
are heavily damped by the Barnes damping process within 
weakly turbulent regimes (Sridhar & Goldreich 1994) However 
the damping weakens in strong turbulence, but according to 
Goldreich & Sridhar | ( 1995), the wave generation of pseudo- 



Alfvenic wavemodes is only possible via three- wave interactions 
by two shear wavemodes. Barnes damping is important for high- 
P plasmas. Since this is not fulfilled for the solar corona, the role 
of pseudo waves should not be ignored. 

Because the model consists only of these two wave types, 
it is suitable to use a description with Alfvenic waves moving 
either forwards or backwards. This is achieved by introducing 
the Elsasser variables (Elsass e7|l950| ) 



w 



v + b- v A e\\ 
v - b + v A e h 

and transforming the Eqs. |5]) into a suitable form of 

i kgkpky 
2~ 



(6) 



(5, - v A k z ) w a = 



k 2 



(d, + v A k z ) w + a = 



2 

2 A- V'fi"i 



i kgk^ky 



wZw^ + w a w 



P' y) 



ah~ + 



(7) 



Obviously, the nonlinearities of Eqs. (|7]) that describe the 
turbulent behaviour of the MHD-plasma cannot be solved in 
Fourier space. Hence the main numerical load is the transfo- 
mation between real- and wavenumber space for each iterative 
step. For this purpose we used the P3DFFT algorithm, which is 



a MPI-parallelised fast Fourier transfomation based on FFTW3 
( |Pekurovsky|20111i| 

One basic problem of spectral methods that use discrete 
Fourier transformation is the aliasing effect. Because of discrete 
sampling in the wavenumber space, high fe-values exhibit errors 
that depend on the structure of the real space fields. Therefore we 
used zero padding, which is also referred to as Orszags 2/3 rule, 
i.e. 2/3 of the wavenumbers below the Nyquist frequency have to 
be truncated to achieve maximum anti-aliasing, hence reducing 
the Fourier space resolution to 1/3 of the original wavenumber 
range (Orszag (1971)). This process is repeated each step, im- 
mediately before calculating the nonlinearities and, accordingly, 
calculating the RHS of the MHD equations. Consequently, the 
change in the antialiasing-range of one MHD-step is physically 
correct, but not the long-term evolution. 

Gismo is capable of using different foward-in-time schemes, 
namely Euler and Runge-Kutta second as well as fourth order. 
All the simulations in this paper has been peformed by using 



Convection 

Dissipation 




Log(k) 

Fig. 1. Possible mechanisms that might influence the evolution 
of an amplified wavemode embedded in a turbulent plasma with 
a Kolmogorov-type power-law energy spectrum. 



To mimic an SEP-event, understanding the underlying mech- 
anisms is crucial. Before we simulated particle scattering at 
modified field-fluctuations, we focussed on the evolution of 
amplified wavemodes within the background turbulence. Since 
streaming particles ejected by the Sun, e.g., protons from coro- 
nal mass ejections, will not sharply amplify a discrete mode 
but a broader interval, a Gaussian distributed shape was as- 
sumed. The streaming instability has the highest growth rate in 
the background field direction. Consequently, we assumed that 
only purely parallel-propagating Alfven waves are modified. 

The Alfven wave generation mechanism by SEPs is not be 
investigated in detail here. The driving mechanism assumed for 
our simulations is the streaming instability. The estimatation of 
the wave growth rate is described in Vainio| ( 2003 ) The stream- 
ing instability is caused by energetic proton scattering off the 
interplanetary Alfven waves. During the the scattering process 
the particle changes its pitch angle cosine by A/i while its mo- 
mentum in the wave frame remains constant. Thus the particles' 
energy in the plasma frame is changed by v A pAfj., consequently 
also happens to the Alfven wave energy due to energy conser- 
vation. Another important instability in the solar corona is the 
electrostatic instability. This is caused by an electron current as 
well as by streaming ions. Ion acoustic waves would be gener- 
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ated by this process. However, for the growth rate of these modes 
a sufficiently high ratio T e /Tf » 1 of the electron and ion tem- 
peratures is crucial. Observations and simulations in the vicinity 
of three solar radii indicate temperature ratios of the order of 
unity (Landi 2 007[|Jin et al.|2 012) In this regime the ion acous- 
tic waves are also efficiently suppressed by Landau damping. For 
further reading about the streaming instability we refer to |Gary| 
(1993), 

These parallel peaked modes are influenced by dissipation, 
diffusion, and convection. An illustration of this is shown in Fig. 
[T] Regarding the evolution within the Fourier space, dissipation 
will damp the wavemode, hence lower the maximum without 
altering the position or broadness of the peak. The dissipation 
of wavemodes is caused by the spatial diffusion term in Eqs. 
(j7]). Convection in Fourier space would shift the position of the 
peak. If diffusive transport in Fourier space were the dominant 
mechanism, it would result in a broader energy distribution. The 
dynamics of convection and diffusion lie within the nonlinear 
terms, therefore one cannot distinguish exactly between the re- 
sponsible terms. 

To investigate the effects of spatial diffusion on the peaks, 
we solved the dissipation equation in wavenumber space 



de 



= -k 2 De 



and calculated its dissipation coefficient 



D = 



1 



k 2 At 



(8) 



(9) 



When the spatial diffusion is the only dissipative effect in 
wavenumber space, this wavenumber dissipation coefficient 
would be the spatial diffusion coefficient. We emphasise that 
diffusion in the wavenumber space is a different process and is 
hence explicitly distinguished from spatial diffusion. In the dis- 
cussed context above the spatial diffusion is clearly connected 
to the wavenumber dissipation. We concentrated our investiga- 
tion on those wavemodes whose peak energy is initiated only. 
For other modes this approach is not feasible because the spatial 
diffusion is not necessarily the dominant process on an arbitrary 
wavemode. 

3. Simulation setup 

To simulate the turbulent plasma in which the SEP-event is set, 
we performed the following type of magnetohydrodynamic tur- 
bulence. 

We used an anisotropic ally driven turbulence with a driv- 
ing range in k-space up to the first five numerical wavenumbers 
in perpendicular (k' ± = 27r[0 • • • 4]) and 15 in parallel direction 
(fcjj = 27r[0 • • • 14]). A remark to the notation, the wavenumber 
is in general defined as k — (2nri)jL where n stands for the 
numerical grid position. For simplicity we used the normalised 
wavenumbers k' = (2n ■ n) throughout. The anisotropy was cho- 
sen for two reasons. First, to mimic the preferential direction of 
the solar wind, where particles that radially stream away from 
the Sun form the Parker spiral. Consequently, these particles can 
deposit their energy in a parallel direction on different scales. 
This is mainly valid in the vicinity of the Sun, in which we are 
interested. Second, a slab-component of SW turbulence is ob- 
served also at small scales in parallel direction. To ensure tur- 
bulence evolution up to high parallel wavenumbers, the driving 
range was extended along the parallel axis. This is necessary be- 
cause the parallel evolution is much weaker than the perpendicu- 
lar. Even though this is primarily a technical aspect to ensure the 



extent of the spectrum to higher k\\, it is still in line with observa- 
tions. An isotropic driver would not yield sufficiently turbulent 
modes at high k\\. The turbulence driving is performed by allo- 
cating an amplitude with a phase to the Elsasser-fields within 
the Fourier space. The amplitude follows a power-law of \k\~ 
and is initialised using a random normal distribution. The phase 
was randomly chosen between zero and 2n. These settings are 
divergence-free and hermitian symmetric. After this initialisa- 
tion the values were scaled to the desired scenario, which in our 
case is a cLB/Bq ratio of roughly 1CT 2 . Note that both species, 
pseudo- and shear Alfven waves, are excited by this type of tur- 
bulence driving, but as presented by Maron & Goldreich (2001) 
the pseudo-waves evolution is strongly suppressed. In this ini- 
tial range energy is injected at discrete times, which leads to a 
saturated turbulence - an equilibrium between dissipation and 
injection. The spatial resolution is 256 3 gridpoints, resulting in 
128 3 points in k-space of which \k'\ = 2n ■ 42 wavemodes are 
active modes that remain unaffected by (anti)aliasing. The hy- 
perdiffusivity coefficient was set to h — 2. 

The simulations of the turbulent background plasma were 
performed assuming an outer scale of L sca i e = 3.4 ■ 10 8 cm. This 
value was estimated using the growth rate from Vainio ( 2003) > 



T(k) = 



2n„v A 



I 



d 3 pvn\k\6{\k\-—)f P , (10) 



with the proton cyclotron frequency a> cp , the proton speed v p , 
the Lorentz factor y, the proton number density n p , fi the pitch 
angle cosine, and f„ the proton distribution function. Here the 
resonance condition for the nth order of interaction 



k v \\ 



nQ. = 0, n e Z 



(11) 



(cf. |Schlickeiser|1989| l was used, where u is the wave frequency, 
and k\\ its parallel wavenumber. Q. is the particle's gyrofrequency 
and V|| its parallel velocity component. Note that Eq. 10 is only 
valid for purely parallel waves and n — +1. Orders of \n]> 1 can 
only be generated by oblique waves. The perpendicular com- 
ponents of the wave would then modify the scattering process 
by nonvanishing Bessel functions. This is discussed in detail by 
Schlickeiser (2002) We used peaks at k — 2n ■ 8 and k-2n- 24, 
which represent proton energies of £ a; 64 MeV and E w 7 MeV 
respectively. Using the resonance condition, this leads to a length 
scale of 



2nn a 
Ascak = —prym p cv » 10 cm. 
eB 



(12) 



The Alfven speed was assumed to be v A = 1.2 ■ 10 8 cm s , 
which leads along with a particle number density of 10 5 cirT 3 to 
a background magnetic field Z?o = 0.174 G. These values resem- 
ble the solar wind environment at three solar radii ( |Vainio et al.| 
(2003)), where particle acceleration by CME-driven shocks is 
strongest. 

The discretisation of the timestep is stable for values up to 
At' — 1 • 10~ n in numerical units or At' = 3.4 • 10~ 3 s for the 
background turbulence. 

If the background plasma simulation reaches the saturated 
state, a Gaussian-distributed energy peak with purely parallel 
k = ke\\ is injected. We chose two different positions of the 
peak in wavenumber space. To investigate the physics of an 
SEP-event, a wavenumber of k\\ = 1.5 ■ 10~ 7 cirT 1 was used. 
This corresponds to a numerical wavenumber of 8, which is 
still within the driving range of the turbulence. The injection at 
smaller scales was represented by a peak at k\\ = 4.4 ■ 10~ 7 cm -1 . 
This value lies at the numerical position 24, which is roughly 
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Table 1. Parameter setup for the simulations. 





^scale 


[cm] 


«rf[cm 3 ] 


B [G] 


v A [cm s _1 ] 


v[num] 


fc-grid 


SI 


3.4 


10 s 


10 5 


0.174 


1.2 ■ 10 s 


1 


128 3 


SII 


3.4 


10 8 


10 5 


0.174 


1.2 - 10 s 


10 


128 3 


SIII 


3.4 


10 8 


10 5 


1.74 


1.2 - 10 9 


10 


128 3 


SIV 


3.4 


10 s 


10 5 


0.00174 


1.2 ■ 10 6 


1 


128 3 



Notes. The wavenumber grid is defined as k = 2^/L sca i c [grid position]. The number density nj connects the background field So with the Alfven 
speed va via ^lAnmrij. 



between the maximum driven wavemode and the anti-aliasing 
truncation edge (which was at 43). We injected the SEP-energy 
gradually over a certain time interval to develop a realistic sce- 
nario. Multiple situations were explored, by using simulations 
with peaks at either position, with either large (growth rate T{) 
or small (growth rate Tz) total amplitude of the Gaussian at the 
final driving step. Because the velocities increase near the peaks, 
the discretisation of the timestep had to be decreased to values of 
At' = 5 ■ 1(T 12 in numerical units or At' = 1.7 ■ 10~ 3 s to sustain 
stability. 

To allow even more diverse case studies, four different ini- 
tial conditions were used, as described in table [TJ In each setup 
a complete evolution of the background turbulence was sim- 
ulated. The first setup SI uses the standard parameters as de- 
scribed above. The main parameters of interest are the resistivity 
of the plasma and the background magnetic field. The results 
in changing these will reveal the effects on the mechanisms de- 
scribed in Fig. [T] An increased resistivity compared to simtype 
SI is approached in SII. A higher value for v is expected to make 
a difference in the spatial diffusion behaviour and would make 
wavenumber dissipation more dominant to the other transport. 
As indicated in Fig. [T] this would lead to a significant damping 
of the peak. The dissipation range of the background turbulence 
will likely be increased by this parameter as well. The third setup 
SIII has a magnetic field increased by a factor of 10. This is to 
examine the influence of a more anisotropic turbulence because 
the perpendicular evolution should be much stronger according 
to Goldreich & Sridhar ( 1995 ) In general, these values may only 
be achieved in magnetic clouds, but it gives valuable information 
on the mechanisms of turbulent transport. The high resistivity is 
necessary because of stability problems with the accompanying 
high Alfven speeds. The last variation of the scenario, which was 
used in SI, is a strongly decreased magnetic background field Bq, 
see simtype SIV. The aim of this artificial scenario is to investi- 
gate strong turbulence at ( » 1 . 



4. Results 

The evolution of the anisotropic background turbulence was sim- 
ulated up to 30 Alfven wave crossing times, which corresponds 
to a simulation of 85 s in physical time. At this point, the tur- 
bulence has reached a saturated state and a Kolmogorov-type 
power-law has evolved over a wide range of wavenumbers (see 
Fig.[2]l. According to |Sridhar & Golckeicl^ ([T994) and |Goldreich 
|& Sri dhar ( 1995}, the 5/3-spectrum is dominant for perpendic- 
ular k ± , whereas the parallel evolution is significantly slower. 
Our simulations confirm this behaviour very clearly. Note that 
the total spectra in Fig. [2] deviate from the Kolmogorov shape, 
especially at high wavenumbers, because the shape of the par- 
allel spectrum is not 5/3. As expected, the spectra are very sen- 
sible to v. A factor of ten increases the dissipation range dras- 
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Fig. 2. Magnetic energy spectrum of the simulated background 
turbulence setups. The plot was made by total integration within 
the Fourier space. 



tically. This is also due to the hyperdiffusivity we used where 
higher wavenumbers are damped by higher power of k (dissipa- 
tion term oc k 4 ) (see Sect. [2}. The dB/Bo ratio of the developed 
turbulence is about 10~ 2 . The magnetic field fluctuations are de- 
fined in Fourier space by 



(dB) 2 



J <?k l -{w-{k) 



w + {k)f. 



(13) 



|*[>0 



The influence of the turbulence strength on the energy transport 
is an important aspect for the peak simulations. We show this 
below. 

Once the background plasma was simulated, a peak is driven 
over a time period of ca. 1.7 s. An exemplary time evolution of 
the peak at normalised wavenumber 24 is shown in Figs. [3] and 
[4] This is a one-dimensional spectrum of the magnetic field en- 
ergy in numerical units that shows a cut along the parallel axis. 
The starting time ? start corresponds to the end of the driving in- 
terval, hence the maximum amplification of the wavemode. The 
timesteps are shown in table[2] For subsequent use the times t^id 
and f en d are introduced. Note that the decay time interval of the 
excited mode at ki = 27T-8 is significantly longer compared to the 
Id — 2n ■ 24 peak. This is again because of the hyperdiffusivity, 
which damps higher modes more strongly because the dissipa- 
tion term is oc k A . The time f en( j is defined at the state where the 
peak has lost nearly its complete energy within the parallel di- 
rection. We took this as the final point of the energy diffusion or 
dissipation. 
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Table 2. Evolution timesteps of the peaks. 





t[s] 


^=2^-8 


^ = 2^24 




h 


= hturt 










h 


— fmid 


51(25.5) 


8.5 


SI 


h 


— ^end 


102 (51) 


17 




h 


= ^start 










h 


— ^mid 


20.4 


1.28 


SII 


h 


— ^end 


40.8 


2.55 




h 


= ^start 










h 


— fmid 


22.53(21.42) 


2.04 


SIII 


h 


— ^end 


45.05 (42.84) 


4.08 




t\ 


— ^start 










'2 


— fmid 


54 (28.9) 


12.75 


SIV 


h 


— ^end 


108 (57.8) 


25.5 





Notes. The labels start, mid and end stand for the times of the decay un- 
til the final dissipation of the peak mode. The time t mK \ is defined as the 
half-time of the decay cycle. Note that the peak at smaller wavenumber 
k'^ = 2n ■ 8 remains visible significantly longer than the other peak. The 
values denoted in brackets are the middle and end of the simulations 
that were not performed until the final dissipation of the peak due to the 
long computational times. In these cases we estimated of the total decay 
by an exponential fit of the decay curve. 



m 1e+10 
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Gaussian fit 
Gaussian fit 
Gaussian fit 
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Fig. 4. Simulation setup SI: One-dimensional energyspectrum 
E(k\{) of the peak at numerical wavenumber 24 (=4.4- 10~ 7 cm -1 ) 
with higher growth rate I^. The influence of the diffusion pro- 
cess is significant because the peak is broadening during time 
evolution. Adjoining maxima develop, e.g. at = 2n ■ 16 and 
38, highlighted by red circles. 
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Fig. 3. Simulation setup SI: Time evolution of a Gaussian- 
distributed amplification at numerical wavenumber 24 (=4.4 • 
10~ 7 cirT 1 ) at the lower growth rate T\, The spectrum is a one- 
dimensional cut along the parallel wavenumber axis where the 
peak is located. The peak is clearly shifted towards smaller k\\. 



Both peaks at Id -2n -24 show broadening of its shape and 
shifting from the initial value to Id = 2n ■ 23 within 17 s, while 
peaks at ki = 27r ■ 8 are only slightly shifted to Id - 2n ■ 7.7 
s within 17 s, but the broadening is clearly visible, as shown in 
Fig.0 

The next step is to investigate the development of the ampli- 
fied wavemodes. If the evolution of the peak were solely gov- 
erned by spatial diffusion, the dissipation coefficient D of Eq. 
([9]) would stay constant in time. Therefore the change in peak 
energy was measured. Two intervals were used for the calcula- 
tion of D, denoted as time-interval T\ and Ti. These intervals are 
different for both peak postions because of the faster decay at 
high wavenumbers. At Id -2n -8 the time intervals are t\ =8.5 
s and T2 = 17 s, whereas at Id - 2n ■ 24 the intervals t\ - 1.7 



it 1e+08 



E3= 18.7 s) 
Gaussian fit 
Gaussian fit 
Gaussian fit 



10 12 14 
k„L/2ji 



Fig. 5. Simulation setup SI: Time evolution of a Gaussian- 
distributed amplification at numerical wavenumber 8 (=1.5 ■ 
1CT 7 cirT 1 ) at the lower growth rate T\. The spectrum is a one- 
dimensional cut along the parallel wavenumber axis where the 
peak is located. Only a slight shift towards small k\\ is observed. 
The broadening is clearly visible, especially on the flanks of the 
Gaussian curve. 



s and T2 = 3.4 s are used. The results of the diffusion coeffients 
are given in table [3] 



To relate the growth rates Fj/2 with the total resulting ampli- 
tude of the Gaussians of the simulations SI-SIII, the energy was 
measured and compared to the background at timestep f stmt each. 
The results are presented in table |4] 

To investigate the direction of the peak evolution in the paral- 
lel and perpendicular directions, two-dimensional contour plots 
were produced. Fig.|6]shows the time evolution of a peak at nor- 
malised Id, = 2ji ■ 8. The two-dimensional spectrum is a contour 
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Table 3. Dissipation coefficients according to Eq. (j^. 





p 


""Peak 


= 2^-8 

1 2 1 


*Peak ~~ 
1 


27T-24 
p 








SlMTYPE SI 






T\ 


5.09- 


10 14 


1.76 • 10 16 4.12 


10 15 


7.92 • 10 15 


T2 


7.23 • 


10 15 


8.79 -10 16 1.55 ■ 


10 16 


2.92 ■ 10 16 








SlMTYPE SII 






T\ 


1.41 ■ 


10 15 


8.33 • 10 15 9.79 • 


10 24 


1.23 • 10 25 


T~2 


3.28- 


10 17 


7.60 • 10 15 4.85 • 


10 30 


4.91 • 10 32 








SlMTYPE SIV 






T] 


1.25 • 


10 14 


6.74 -10 13 5.84- 


10 14 




T2 


4.43 • 


10 14 


8.00 • 10 13 1.47- 


10 16 





Notes. The coefficients are calculated for different simulation types (SI, 
SII, SIV) at different time intervals (ji,T?), for different growth rates 
(ri,T2), and for the two driving wave numbers (k = 2n ■ 8, 27r • 24). 
Not calculated are the coefficients for SIII because two parameters were 
changed in this setup. Therefore it is not directly comparable to the other 
simulations. The diffusion coefficients are given in cm 2 s~'. 



Table 4. Energy deposited at the driven peak wave number. 



r 


^T-eak 


= 2n- 8 

r 2 


r 


k' 

""Peak 


= 2n ■ 24 

r 2 








SlMTYPE SI 






1.90 


10 7 


1.95 • 


10 9 1.66 


10 18 


2.29 • 10 20 








SlMTYPE SII 






9.17 


10 7 


9.22- 


10' 2.10- 


10 18 


2.10- 10 20 








SlMTYPE SIII 






4.41 


10 7 


4.39- 


10" 4.98 • 


10 19 


4.98 ■ 10 21 








SlMTYPE SIV 






1.45 


10 6 


1.54- 


10 8 5.17 


10 9 





Notes. The energy is given as the ratio of is^Peak, fmax)/£(£peak, fstan)- 



plot of the power spectral density of the magnetic field that was 
calculated by cylindrical integration in fc-space. The single con- 
tours are scaled logarithmically. 

During the evolution, higher harmonics of the initial peak 
arise. To investigate these in greater detail, we measured the en- 
ergy of the initial peak and its first harmonic. The result is shown 
in Fig. [7] Note that the generation of these modes starts at higher 
perpendicular wavenumbers (see Fig. |6j most picture left, first 
harmonic at k ± m 15) and not at purely parallel k. 

In addition to the observed higher harmonics, another basic 
result of the simulation type SI is the discrete generation of other 
modes along the k\\ axis especially at higher amplitudes. As seen 
in Fig. [4] adjoining maxima develop next to the main Gaussian 
curve, e.g. at positions Id = 2n ■ 16 and 38. The way the peak 
develops appears to vary with the amplitude, the peak with Ti 
generates clearly fewer other modes than the larger peak with Yj- 
Both peaks show a significant change of thei original wavenum- 
ber position and a clear broadening. The two-dimensional spec- 



Energy of the 1st harmonic 
Fit x 2 



1e+12 1e+13 1e+14 1e+15 1e+16 1e+17 

E initial peak( k ||) [numerical units] 

Fig. 7. Energy dependency of the first harmonic on the initial 
peak energy. A fit resulted in a quadratic function. The simula- 
tion setup SI was used. 



tra reveal a strong perpendicular development in k-space espe- 
cially at large \k\. During the decay the evolution becomes a lit- 
tle more isotropic, but the preferential direction of evolution is 
clearly perpendicular. 




Fig. 8. Simulation setup SI: Comparison between growth rates 
Fi and T2 at peak position k',, = 2n ■ 24 at maximum drive time. 
Although the effective energy input varies by a factor 100 (see 
table [4} another transport mechanism becomes dominant. The 
smaller peak (Fi) develops dominantly in the perpendicular di- 
rection while the evolution of bigger peak is more isotropic 
and tends towards smaller k\\. The colours indicate the logarithm 
of the total spectral energy. 



The importance of the amplitudes T1/2 is also visible in Fig. 
[8] The development of each peak is very different. An interesting 
result is a more dominant evolution of the Id - 2n ■ 24 peak 
with the high growth rate F2 that is towards smaller k\\ directed. 
In the peak-dominated region the turbulence seems to increase. 
Therefore the f parameter (Eq. |2]i is of interest. Fig. [9] shows a 
map of values of £ = [0.01 • ■ ■ 0.15], i.e. near the critical balance. 
The same plot for the lower growth rate would be empty. This 
also indicates that the high ( values along the A; ± -axis in Fig. [9] 
stem from interactions with the peak. 

The peaks at Id — 2n ■ 8 remain much longer than the higher 
modes, see table [2] This is due to the higher k on which the 
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10 20 30 40 10 20 30 40 10 20 30 40 



k H k H k t 

Fig. 6. Two-dimensional magnetic energy spectra of a peak at normalised wavenumber 8. Red regions are at higher energies com- 
pared to the blue ones. The parallel and perpendicular wavenumbers are given as absolute values. The time development is shown 
for mid-drive (Af w 0.85 s), max-peak (Af * 1.7 s) and the decay 17 s after the driving. The colours of the contours were normalised 
for comparison between the three plots. The colours indicate the logarithm of the total spectral energy. The simulation setup SI was 
used. 





10 20 30 40 




20 30 



Fig. 10. Two-dimensional peak evolution in setup SII at maxi- 
mum drive time. The higher harmonics develop at lower k ± com- 
pared to SI. The higher growth rate (T2, right panel) shows a 
strong parallel evolution. The smaller amplitude (T\, left panel) 
develops dominantly towards higher k ± . The colours indicate the 
logarithm of the total spectral energy. 



Fig. 9. Map of the critical balance parameter f for the right-hand 
plot in Fig. [8] The contours linearly represent values between 
0.01 and 0.14 of the integral values of £(k\\,k ± ). The peak struc- 
ture and values near the A: ± -axis are clearly visible. 



dissipation process depends. Furthermore, the hyperdiffusivity 
damps higher modes more strongly (oc k A ). 

The simulation of the same peaks in setup SII with higher 
v reveals one key feature. The shift towards smaller k\\ positions 
occurs on both, k',. = 2n-% and 24 peaks and is stronger compared 
to SI. The shifting is significant, especially at k^ = 2n ■ 24, which 
is shown in Fig[TT] Within 2.55 s the original position changes 
from Id = 2n ■ 24 to 22. The peak amplitude is again impor- 
tant for the evolution. We observed that higher growth rates lead 



to an isotropic evolution, whereas the lower growth rates show 
strongly perpendicular development. The effective peak energy 
is lower and consequently the energy transport towards higher 
wavenumbers is more restricted. There are also fewer of higher 
harmonics. This can be observed by direct comparison of the left 
plot in Fig. 10 with the middle plot in Fig. [6] As expected, the 
decay of the energy is faster compared to SI (see table [2]). 

The simulation setup SIII reveals very strong perpendicular 
evolution for all peaks. Two examples are given in Fig. [12] The 
growth rates seem not to have a strong influence on the develop- 
ment. Only more higher harmonics of Id = 2n- 8 peak are visible 
with r%. The time of total energy decay of the peaks is slightly 
longer than in the simulations SII. 

In the last simulation setup SIV the magnetic background 
field was reduced by two orders of magnitude. The resulting 
dB/Bo ratio is of the order of 10 and the turbulence develop- 
ment is highly isotropic. The peaks at Id = 2ji ■ 8 and k!, -2n-2A 
both at growth rate F2 show very interesting features. 
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Fig. 11. Simulation setup SII: Time evolution of a Gaussian dis- 
tributed amplification at ki — 2n ■ 24 at the smaller growth rate 
T\. Again a ID cut spectrum along k\\ axis is shown. The shift of 
the peak is stronger compared to SI. 



■i 20 





Fig. 12. Evolution of the two peak positions at f m id in setup SIII. 
Very strong perpendicular development in all SIII simulations is 
observed. The edge of the parallel driving range at Id — 2n ■ 14 
is clearly visible in the right-hand panel. 




Fig. 13. Two-dimensional peak evolution in setup SIV. Left: 
Development of the peak at ki = 27r ■ 8 with growth rate L^. 
The figure presents the state 1.28 s after the maximum driven 
peak (to). Right panel: Corresponding map of the critical bal- 
ance parameter f . Each contour represents integral values above 
^ = 0.1. The colour scaling is linear. 




30 40 




30 40 



Fig. 14. Two-dimensional peak evolution in setup SIV. Left: 
Development of the peak at k^ = 2n ■ 24 with growth rate F\. 
The figure presents the state 0.68 s after the maximum driven 
peak (fo). Right panel: Corresponding map of the critical bal- 
ance parameter Each contour represents integral values above 
£ = 0.1. The colour scaling is linear. 



We discuss possible explanations for these phenomena be- 



In addition to the typical peak evolution and generation of 
higher harmonics, other structures also arise at high k x . As 
shown in the Figs. 13 and 14 both peak positions generate 
these structures, but at various positions. During the develop- 
ment of the peak at position Id = 27r ■ 8 theses structures arise 
perdominantly at high (k'^k'j) locations, e.g. (2n ■ 10, 2tt • 18), 
(2n ■ 10, 2n ■ 38), (2tt • 18, In ■ 32) and (2n ■ 24, 2n ■ 18), see Fig 



13 The structures evolving with the ki = 2n ■ 24 simulations 
are instead located at middle k ± but low k\\, e.g. (27r • 5, 2n • 13), 
(2n ■ 5, 2n ■ 17) and (27r • 8, 27r • 22), see Fig 14 In contrast to 



the higher harmonics, these structures are not necessarily inte- 
gral multiples of the initial peaked mode. A map of the critical 
balance parameter f was calculated for both plots. Around the 
wavenumber of the peaks and along the k ± axis this parameter 
is of order 1. Interestingly, during the k'« - 2n ■ 24 simulation 
this parameter increases, especially at the k ± axis. In contrast to 
the other setups, the higher harmonics are located along the k\\ 
axis at low or zero k x , A significant shift along this axis towards 
smaller parallel wavenumbers is observed, e.g. as shown in Fig. 

m 



low. 



5. Discussion 

As discussed in Sect. 2] there are three possibilities how an ex- 
cited wavemode can develop: through diffusion, dissipation and 
convection. 

A general conclusion from our simulations is that the dy- 
namics of these mechanisms are strong at high wavenumbers. 
Especially for the dissipation process this is not unexpected be- 
cause of it strongly depends on th wavenumber. The hyperdiffu- 
sivity might amplify this effect because of the higher power in k. 
The Figs. [3] and |4] clearly show a rapid dissipation of energy be- 
cause the peak loses amplitude very fast. Also table ^indicates 
a decay in short timescales at high wavenumbers. However, a 
broadening also arises at the flanks of the Gaussian distribution. 
The broadening is strong for the higher growth rate F2. Within 
a time interval of 1 .7 s after the driving range the FWHM is in- 
creased by roughly 70% at F2, whereas the peak with smaller 
amplitude is broadened by ca. 15 %. A possible explanation is 
the equilibrium between energy- and enstrophy cascade, which 
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causes a similar flow of energy to large and small wavenumbers 
( |Mininni & Pouquet|2009] l 



si, t-w, 

SI, t.17s 
SI", t-tsa, 
Sill. t.17s 

siv, l-W, 

SIV. t.17s 



i ! 



. VV 



Fig. 15. Time evolution of the enstrophy for the simulation se- 
tups SI, SIII and SIV. The strong magnetic background field of 
SIII prevents the enstrophy from developing while a clear change 
in SIV is visible. 



The convective mechanism shifts the maximum of the peak. 
Convection is a slow process compared to dissipation and dif- 
fusion within the simulated regime. Nevertheless we were able 
to observe it within our simulations, e.g. in Fig. [3] The transport 
is towards smaller wavenumbers, which indicates an enstrophy 
cascade. This effect is more typical for two-dimensional plas- 
mas. The MHD-development in anisotropic plasmas is mostly 
effectively two-dimensional. This lead to inverse energy cas- 
cades as well as upscaling enstrophy cascades. These mecha- 
nisms generate larger vortices and consequently transfer energy 
to smaller wavenumbers. 



SI: E_„(t= 17 s) 
SIII: E™ t= 17 s 
SIV: E3-17sj 
Gaussian fit 
Gaussian fit 
Gaussian fit 
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Fig. 16. Comparison of the convective peak shifting between the 
setups with changing background field Bo SI, SIII and SIV. The 
strong Bo in SIII causes no observable shifting at all, whereas 
the weak Bo leads to a significant change of the original peak 
position. One possible explanation is an enstrophy cascade. 



To investigate this behaviour in more detail, we concentrated 
on the peak Id — 2n ■ 8 with growth rate H within the setups 
of changing magnetic background field SI, SIII and SIV. The 
comparsion of these setups is shown in Fig. 16 We observe a 



slight shift of the in setup SI after 17s of the maximum drive 
from position 8 to 7.7. The simulation SIII with increased Bq 
shows no transport of the peak. After the same time interval it 
is still at position 8. The evolution in setup SIV with a small 
Bo is very strong. After 17s the peak has shifted by roughly 1.5 
grid positions from k^ — 2n ■ 8 to 6.5. The development of an 
enstrophy cascade explains this behaviour. The strong magnetic 
field effectively makes the MHD-evolution one-dimensional. As 
Bo decreases, the evolution becomes less restricted to the mag- 
netic background field. Consequently, the enstrophy cascade in- 
creases. Fig. [15] clearly supports this explanation. The enstrophy 
was calculated by 



s(k) 



-I 



<?k \kxu\ 2 . 



(14) 



All two-dimensional spectra show a strong evolution in the 
perpendicular direction. This is consistent with the theories of 
ISridhar & Goldreich| fl994) and |Goldreich~&"Sr idhar ( 19951 
within a turbulent plasma, as explained in Sect. [2] The process 
of this perpendicular evolution is caused by the energy cascade. 
Mainly Fig. 12 shows strong perpendicular behaviour, whereas 
the evolution clearly becomes more isotropic in SIV (see Figs. 
[13] and [14 1. This is due to the increasing strength of the turbu- 
lence for the cascade, which is expressed by the dB/Bo ratio. 

The dissipation coefficients presented in table [3] do remain 
constant especially for all simulations with F2 at Id = 27r • 8. This 
implies that spatial diffusion is the dominant process for these 
simulations. The most significant change of the dissipation coef- 
ficient is at Id — 2n-24 for setup SII, between T\ and t^- We inter- 
pret this as a nonlinear effect where wavemodes are generated, 
which in turn triggers the cascade. This influences the energy of 
the initialised Gaussian very strongly. As expected, the dissipa- 
tion coefficients for setup SII are stronger for the k^ = 2n ■ 24 
peaks. The spatial diffusion coefficient is connected with the 
wavenumber dissipation process via v. This is because the spa- 
tial diffusion process is the only mechanism that leads to energy 
losses in the k space (see Sect.[2|. This is at least valid for wave- 
modes far below the antialiasing edge, where energy is artifi- 
cially removed as well. A connection of wavenumber diffusion 
to v is observed in Fig. 10 where the Gaussian shape is signifi- 



cantly broadened, which is caused by the diffusion process. 

The observed higher harmonics could resemble a three-wave 
process k% + k% — > fci6. This is supported by Fig. [7] The en- 
ergy dependency between peak and the next higher harmonic 
is quadratic. This is also the case for three-wave interactions 
(Sagdeev & Galeev 1969) As pointed out before, this is for- 
bidden for Alfven waves by wave interaction processes. Just 
oppositely directed wavepackages can collide, hence momen- 
tum conservation would be violated by the proccess described 
above. This wave interaction can only take place with oppo- 
sitely directed waves of the background plasma. This must also 
be true since a cross-check simulation of the peak without turbu- 
lent background did not show these harmonics. Another expla- 
nation is given by Galinsky et al. ( 19971): Alfven waves interact 
with themselves, which leads to wave steepening. 

Investigating the strong turbulence evolution in simulation 
SIV shows unexpected structures at high k x (see Fig. 13 and 



14 1. This effect might be caused by the critical balance of strong 
turbulence within the Goldreich and Sridhar description. This 
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requires the parameter £ to be ~ 1, which is the case for loca- 
tions in the vicinity of the peaks and along the perpendicular 
axis. Mainly the Id = 2n ■ 24 peak seems to amplify the re- 
gion at k'j_ — 2n ■ 12 and Id = 0. Values of £ > 0.1 lead to the 
development of these structures, but it is not possible to con- 
clude whether the stuctures arise first and then generate higher 
values of or vice versa. Nevertheless, the turbulence strength 



increases within these regions, which agrees with |Goldreich & 
Sridhar ( 1995 ) where ( is assumed to become unity during the 



turbulence evolution. In addition to setup SIV, SI also shows this 
behaviour as presented in Fig. [8] and [9] More investigation is 
needed to clarify the underlying processes. 



6. Conclusion 

We analysed the evolution of waves generated by proton beams 
in a turbulent medium. This evolution may play an important 
role in diffusive shock acceleration in the heliosphere. Our study 
has revealed that three different processes as sketched in Fig. [T] 
are taking place. 

The most interesting question is wheter wavemodes excited 
by particles of a certain energy yield wavemodes that interact 
with other particle energies. The observed shifting of the initial 
parallel wavenumber position towards smaller k\\ influences the 
particle acceleration at these modes. As shown in Eq. 1 1 and in 



the corresponding section, this means that higher energetic parti- 
cles can be accelerated because the wavemode develops towards 
higher spatial scales ( Vainio & Laitinen 2007 ) The shift towards 
smaller k\\ is indeed fairly minor in this simulation. This is be- 
cause of the injection of energy at only one single wavenum- 
ber and the limited simulation time. We also expect an effect 
through nonlinear amplification: Waves at smaller k\\ accelerate 
higher energetic particles, again injecting energy at lower k\\. On 
the other hand a strong evolution towards high k\\ has also been 
observed in terms of development of higher harmonics of the 
initialised mode. Consequently, particle acceleration at lower en- 
ergies is also possible. It should be noted, however, that this is 
not consistent with isotropic diffusion of energy in wavenumber 
space as assumed in Vainio & Laitinenj ( 2001) This means that 
the previous models of the wave-particle system will have to be 
updated accordingly to account for the strong dependence of en- 
ergy transport on the direction in wavenumber space. Especially 
the development of the higher harmonics contradicts an identical 
forward and backward cascade. 

The strong perpendicular evolution of the peak initialised 
with purely parallel propagating waves causes higher orders 
(\m\ > 1, see Eq. \TT\ of resonance between solar particles and 
the amplified mode. This is because Eq. 10 has to be modified in 
this case ( Schlickeiser 2002) 

Owing to limited computational power we have not been 
able to investigate the effect of critical balance in detail. But we 
note that under certain conditions (k\\/k ± , amplitude) the evolu- 
tion may be governed by the critical balance. 
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